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Abstract 

The computational method developed to study the cryo- 
genic fluid characteristics inside a fuel tank in a hyper- 
sonic aircraft is presented. The model simulates a rapid 
draining of the tank by modeling the ullage vapor and the 
cryogenic liquid with a moving interface. A mathematical 
transformation was developed and applied to the Navier- 
Stokes equations to account for the moving interface. The 
formulation of the numerical method is a transient hybrid 
explicit-implicit technique where the pressure term in the 
momentum equations is approximated to first order in time 
by combining the continuity equation with an ideal equa- 
tion of state. 

Nomenclature 

a Van der Waals constant, Pa m 6 

b Van der Waals constant, m 3 

c speed of sound, m/sec 

FI flux, radial direction 

FJ flux, circumferential direction 

g gravity, m/sec 2 

h enthalpy, J/kg 

k conductivity, J/m sec K 

MSF Mach scaling factor, c rea] /c numerical 
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P pressure, Pa 

q heat flux, J/m 2 sec 

R radius relative to interface, f (0, f), m 

R fluid constant, J/kg K 

R ideal fluid constant, J/kg K 

r radial coordinate, m 

r transformed radius coordinate, r/R 

T temperature, K 

t time, sec 

V velocity, m/sec 

A r distance between radial grid points, m 

Ai3 distance between circumferential grid points, 

rad 

a thermal diffusivity, m 2 /sec 

p viscosity, kg/m sec 

v specific volume, m 3 /kg 

0 circumferential coordinate, rad 

p density, kg/m 3 

transformed circumferential coordinate, rad 

Superscripts 

n time step index 

Subscripts 

i, j spacial coordinate indexes 



i +, i spacial coordinate at i + i - 

j +, j spacial coordinate at j + j - 

r, r radial direction 

s constant entropy 

0,d circumferential direction 


1 

2 

1 

2 


Introduction 


Liquid hydrogen has long been considered one of the 
most advantageous fuels for hypersonic aircraft. The 
National Advisory Committee on Aeronautics (NACA), 
Lewis Flight Propulsion Laboratory, Cleveland, Ohio, 
began to study hydrogen as a fuel source in the 1950’s. 
The NACA study included a flight test program with a 
U.S. Air Force B-57 which had one engine modified for 
operation with hydrogen fuel. 1 The test program was suc- 
cessful, and the engine typically operated at an altitude of 
50,000 ft and Mach 0.75 for 17 min with hydrogen fuel. 
Since this initial test program, there have been several 
studies to investigate the aspects of liquid hydrogen- 
fueled aircraft. An excellent review is presented by 
Brewer. 2 

Developing reusable flight-weight cryogenic fuel tanks 
is one of the most difficult challenges in designing 
advanced hypersonic aircraft and the next generation of 
spacecraft. The most recent cryogenic fuel tank studies 
were conducted for the National Aero-Space Plane. 
McDonnell Douglas Corporation, St. Louis, Missouri, 
constructed and tested with liquid hydrogen a 3600-liter, 
graphite-reinforced, epoxy thermoset composite tank. 3 
The General Dynamics Space Systems Division, San 
Diego, California, also constructed and tested a 2600-liter 
graphite and epoxy liquid hydrogen tank. 4 Presently, 
NASA Dryden Flight Research Center, Edwards, Califor- 
nia, is involved in a research program to continue the 
study and to increase understanding of liquid hydrogen 
cryogenic fuel tanks. The initial phase of the program 
includes the design of a test facility and the development 
of a generic research cryogenic tank (fig. 1) for initial test- 
ing and verification of analytical models. Testing of the 
tank will simulate the severe aeroheating effects of hyper- 
sonic flight by subjecting the tank exterior to large heat 
fluxes. Details of the generic research cryogenic tank and 
test setup have been reported. 5 

These hypersonic cryogenic fuel tanks must be able to 
survive the thermal gradients of exterior gaseous aeroheat- 
ing temperatures as high as 7000 K and internal cryogenic 
temperatures as low as 20 K. See reference 6 for a detailed 
discussion on hypersonic aerodynamic heating. The 
understanding of the cryogenic fluid dynamic motion 


inside the tank is essential to predict and analyze the ther- 
mal stress gradients though the tank wall. The fluid 
dynamic model will predict cryogenic fluid temperatures, 
heat-transfer coefficients on the inner tank wall, and strati- 
fication effects in the ullage gas. An analysis using 
SINDA85/FLUINT, which showed that these are the rele- 
vant interior parameters in determining the fuel tank ther- 
mal stresses, has been performed. 5 

The generic research cryogenic tank which this model 
was developed to simulate has met with budgetary con- 
strains, and test data for verification of the model are 
unavailable. Comparison of the model with test data will 
be the subject of future efforts. This paper presents the 
computational fluid dynamic modeling development effort 
to simulate the generic research cryogenic tank when sub- 
jected to high external heat flux and rapid draining. The 
numerical methods were developed with emphasis on 
increasing the time step to reduce computer resources. 

Model Description 

Figure 2 shows the cross-sectional area of the generic 
research cryogenic tank modeled two dimensionally in 
polar coordinates. The interface between the cryogenic liq- 
uid and the ullage vapor is depicted as a dark horizontal 
line. This interface moves down as the tank is drained of 
liquid. The liquid exits the model through the bottom drain 
port, and ullage gas enters the model through the top pres- 
surization port. The polar coordinate system was chosen to 
be relative to the moving interface to facilitate a coordi- 
nate transformation. The model employs the full set of 
Navier-Stokes equations with the exception that energy 
dissipation is neglected in the energy equation. For com- 
pleteness, the applicable equations are given below. 
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Boundary conditions for the model are divided into 
three regions: tank wall, ullage pressurant gas port, and 
liquid drain port. The wall boundary conditions are 
defined as the no-slip condition, and the heat flux is speci- 
fied as a function of circumference and time. These condi- 
tions are expressed as follows: 


v r = v e = 0 


q = /(6,0 

The ullage pressurant gas port boundary is specified at 
constant pressure and temperature with no heat flux: 


P = P 


gas 


T = T 


gas 


q = o 

The liquid drain port boundary is specified as a free 
boundary with no circumferential velocity, at constant 
pressure, and with no heat flux: 

v e = o 


P ~ P liquid 

q = o 

The interface movement is calculated from the liquid drain 
rate and the geometry. Interface conditions are equal tem- 
perature and equal velocity: 

T =7 

gas liquid 

= y e 

Vgas ° liquid 


Because of the numerical difficulties associated with mov- 
ing boundaries, a transformation is defined which allows 
the transformed computational space to be invariant to the 
moving vapor and liquid interface. This transformation is 
graphically shown in figure 3 and is defined as follows: 
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This transformation is now applied to the model equa- 
tions to obtain the applicable Navier-Stokes equations for 
the computational space. These equations are presented as 
follows: 
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Gas and liquid states remain unchanged. 

Boundary conditions at tank wall 

V = V, = 0 

q = mo 

The boundary conditions at the ullage pressurant and 
liquid drain ports remain unchanged. Interface conditions 
also remain unchanged. 

Note that the transformed equations reduce to the stan- 
dard Navier-Stokes equations for cylindrical coordinates 
when 
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Numerical Method 


Finite-difference equations were specifically developed 
for this model using a volume integral method technique 
which is discussed in lecture notes by Lick. 8 The develop- 
ment of the continuity and energy equations is purely 
explicit, time-marching algorithms. The development of 
the momentum equations is longitudinally implicit, trans- 
versely explicit, time-marching algorithms. Longitudinally 
implicit means, for example, that in solving for V r in the 
radial momentum equation all terms that are discretized in 
the radial direction are implicit. Conversely, terms that are 
discretized in the theta direction are explicit (reverse for 
the theta momentum equation). 


The longitudinally implicit approach to the momentum 
equations was chosen to increase the stability of the solu- 
tion to the nonlinear Navier-Stokes equations. This 
approach improved Reynolds cell instability associated 
with purely explicit time-marching techniques. The finite- 
difference equations for the transformed Navier-Stokes 
equations are second-order accurate in space and first- 
order accurate in time and are presented as follows: 
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where the pressure term, P , in the radial momentum 
equation is approximated to first order in time by combin- 
ing the continuity equation with an ideal equation of state. 
The result is as follows: 
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complete set of time-marching, finite-difference equations. 
The first step in the numerical procedure is to solve the 
momentum equations for the velocities at the next time 
step, V ” + and V# + \ The formulation of the momen- 
tum difference equations is such that a tridiagonal matrix 
is solved for each ray of radial velocities. In addition, a 
tridiagonal matrix is solved for each circlet of circumfer- 
ential velocities. The average velocity over the time step, 
V avg = (y n + l + V n )/2 , is then used in the continuity 
and energy equations to obtain density and temperature^ for 
the next time step, p " and T ” j . The pressure, P i j , is 

then calculated from the state equation as a function of 

n + I , — n + 1 

pm ^ r i, j ■ 

A staggered grid was used to develop the circumferen- 
tial and radial velocities (that is, velocity nodes centered 
between pressure nodes). The staggered grid was selected 
to avoid special consideration of the no-slip radial bound- 
ary condition. This special consideration arises because 
longitudinal pressure disturbances or acoustical waves 
traveling in the radial direction are not correctly modeled 
in a standard grid by assuming the wall velocity to be zero. 
The fluid molecules in the vicinity near the wall have 
velocity when disturbed by radial longitudinal acoustical 
waves. Thus, specifying a zero velocity for a wall bound- 
ary node does not represent the total average velocity for 
the boundary node when using a standard grid. This diffi- 
culty is avoided by using a staggered grid and a no-flux 
boundary condition at the wall because the radial veloci- 
ties are not calculated at the boundary in the staggered grid 
but are calculated at a location of Ar/2 from the boundary. 
In addition to the staggered grid, variable node spacing 
was incorporated to model the tank wall boundary layer 
with more densely populated nodes. 
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The boundary and interface finite-difference equations 
are derived directly from the volume integral method tech- 
nique and are very similar to these equations. Because of 
the similarity, they are not shown. These equations form a 


The stiffness often associated with the Navier-Stokes 
equations can present numerical difficulties. Stiffness is a 
term which characterizes a system which has two vastly 
different time scales. In this problem, the stiffness results 
from the disparity between the numerical time step 
associated with the speed of sound in the fluid, of the order 
of 0.00001 sec, and the drain time of the tank, of the order 
of 1000 sec. The numerical solution would require an 
order of 100,000,000 time steps. One method of dealing 
with the difficulties associated with stiffness is to adjust 
the characteristic time scales of the variables. 9 - 10 The stiff- 
ness in this problem is reduced by adjusting the Mach 
number of the fluid by artificially decreasing the speed of 
sound in the computation. The Mach number of the flow 
within the cryogenic tank is expected to be of the order of 
0.0001 which is nearly incompressible. Because com- 
pressible flows with Mach numbers of 0.1 or less can be 
approximated as incompressible flows with less than 
1 percent error 11 , it is reasonable to expect that the Mach 
number of this problem can be increased from 0.0001 to 
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0.1 with less than 1 percent error. The Mach number is 
increased in the computation by scaling the pressure vari- 
able as it relates to the speed of sound through the relation, 



The numerical time step associated with the scaled pres- 
sure variable can be increased to the order of 0.01 sec, 
thus reducing the required time steps for solution to 
100,000. 

Results and Discussion 

These numerical calculations were conducted on small 
grids (0.0 1-m radius with 280 nodes and 0.05-m radius 
with 345 nodes) to save computer resources during model 
development. Modeling of the generic research cryogenic 
tank will require a larger grid (0.75-m radius with approx- 
imately 5,000 to 10,000 nodes). This grid will substan- 
tially increase computer resources. Full-scale calculations 
were considered an unwarranted use of computer 
resources because testing of the generic research cryo- 
genic tank has not begun. Data for model comparison are 
currently not available. 

Figures 4 and 5 show the numerical effect of scaling the 
Mach number. These figures present time histories of 
velocity at the center of the gas region (fig. 4) and at the 
center of the liquid region (fig. 5) for various Mach scaling 
factors (MSF). The calculations were performed assuming 
a tank radius of 0.01 m with a 10- by 28-node grid. Ini- 
tially, the gas and liquid regions are at rest. The model is 
begun by applying a constant impulsive drain rate velocity 
at the gas and liquid interface. These figures represent the 
time response to an impulsively applied drain rate. Excel- 
lent agreement for Mach scaling factors up to 100 and fair 
agreement for Mach scaling factors up to 1000 are also 
shown. The time step required for stability in these calcu- 
lations were l.e-9 sec (MSF =1), l.e-7 sec (MSF = 100), 
and l.e-6 sec (MSF = 1000). Thus by using a Mach scal- 
ing factor of 1000, computational time is decreased by 
three orders of magnitude in this test case. 

Figures 4 and 5 also show the physical modeling effect 
of scaling the Mach numbers. Note that the time axis is a 
ratio of the Mach scaling factor. Thus when an impulsive 
disturbance is applied to the model, scaling the Mach 
number increases the time required for the modeled fluid 
to respond to the disturbance. The final effect of the distur- 
bance, or the velocities, remain unchanged. In this model, 
the interest is in fluid characteristics which occur over 
large time scales (minutes). As a result, characteristic 
responses which occur in time scales of microseconds can 
be increased to seconds as long as they do not affect the 
final result. 


Figures 6 and 7 show the numerical results of applying a 
first-order correction to the pressure term 



in the momentum equations for the ullage gas and liquid 
regions, respectively. In these figures, R/R* is the ratio of 
the fluid constant used in the numerical calculations to the 
ideal fluid constant. Thus, an R/R* = 0 would imply no 
first-order correction. An R/R* = 1 would imply an ideal 
first-order correction. The deviation is defined as the max- 
imum velocity difference between a particular test case 
and a standard test case calculated with a time step of l.e- 
6 sec and R/R* set at zero. 

These calculations were performed assuming a tank 
radius of 0.01 m with a 10- by 28-node grid and a Mach 
scaling factor of 100. As expected, the first-order correc- 
tion term increased the stability of the program and 
allowed larger time steps in the calculations. However, the 
large deviations with increasing R/R* for the liquid (17 
percent at R/R* = 1, Af = 0.5e-5) were unanticipated. The 
large difference in deviation between the gas and liquid 
regions may result from the difference in Reynolds cell 
number. The average Reynolds cell number in the gas 
region is 0.65 compared to 8.5 in the liquid region. Rey- 
nolds cell numbers greater than 2 are known to cause 
instabilities in the finite-difference technique. An R/R * = 
0.25 appears to be the optimal value for this test case. As a 
result, the computational time could be reduced by half 
(by doubling time step) with a loss in numerical accuracy 
of less than 3 percent in the gas region and 10 percent in 
the liquid region. 

Figures 8(a), (b), and (c) show the effect of the moving 
interface. These figures present velocity vectors at two- 
thirds, one-half, and one-third of the gas and liquid inter- 
face heights. The calculations were conducted on one side 
of the tank as symmetry is assumed, a tank radius of 0.05 
m, a 15- by 23-node grid, a Mach scaling factor of 100, an 
R/R* set at 0.25, and a drain rate of 0.002 m 3 /sec. 


7 


Concluding Remarks 

A computational fluid dynamic model has been devel- 
oped to model a cryogenic fluid tank while subjected to 
rapid drain rate and high heat flux. The model incorporates 
a moving interface, a transient hybrid explicit-implicit 
finite-difference technique, a Mach scaling factor, and a 
first-order correction to the pressure term in the momen- 
tum equations. The Mach scaling factor decreases the 
computational time of the model without a significant loss 
of accuracy. The first-order correction to the pressure term 
increases stability and decreases computational time with 
some loss of accuracy. 

The generic research cryogenic tank which this model 
was developed to simulate has met with budgetary con- 
strains, and lest data for verification of the model are 
unavailable. Comparison of the model with test data will 
be the subject of future efforts. 
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Pressurization gas: Specify incoming 



Figure 2. Modeled cross-sectional area. 
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Figure 3. Coordinate transformation relative to the vapor and liquid interface. 


9 



Velocity, 

m/sec 


Velocity, 

m/sec 



Theta 

velocity, 

MSF 

1 

10 

100 

1000 

10000 


♦ 

▲ 


Radial 

velocity, 

MSF 

1 

10 

100 

1000 

10000 


Figure 4. Numerical effect of scaling the Mach number in the gas region. 
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Figure 5. Numerical effect of scaling the Mach number in the liquid region. 
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Figure 6. Numerical deviation of first-order correction to the pressure term in the gas region. (R/R* - 0 implies no first- 
order correction. R/R* = 1 implies an ideal first-order correction.) 



Figure 7. Numerical deviation of first-order correction to the pressure term in the liquid region. 
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(b) One-half. 

Figure 8. Velocity vectors at various liquid levels. 
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